Setup

set.seed(234876)

library(MCMCglmm)
## Loading required package: Matrix
## Loading required package: coda
## Loading required package: ape
library(tidyverse)
## ── Attaching packages ──────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.0.0     ✔ purrr   0.2.5
## ✔ tibble  1.4.2     ✔ dplyr   0.7.6
## ✔ tidyr   0.8.1     ✔ stringr 1.3.1
## ✔ readr   1.1.1     ✔ forcats 0.3.0
## ── Conflicts ─────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
library(parallel)
library(knitr)
# purl("07_h2surv_multivariate.Rmd", output = "07_h2surv_multivariate.R", documentation = 2)

Data and functions for color mapping

source("color_map.R")
source("ggplot_theme.R")
source("h2life_load_data.R")
## Parsed with column specification:
## cols(
##   setDate = col_date(format = ""),
##   flipDate = col_date(format = ""),
##   age = col_integer(),
##   NewAge = col_integer(),
##   fID = col_character(),
##   id = col_character(),
##   sireid = col_character(),
##   damid = col_character(),
##   repl = col_character(),
##   treat = col_character(),
##   NstartF = col_integer(),
##   status = col_integer()
## )
h2life$treat <- str_replace(h2life$treat, "LY", "DR")
h2life$treat <- factor(h2life$treat)
h2life$treat <- relevel(h2life$treat, "STD")
genet_corr <- tibble(model = character(),
                     HS_STD = numeric(),
                     DR_STD = numeric(),
                     HS_DR = numeric(),
                     n_eff = numeric())

iter <- 1e6
burnin <- 2e4
thinning <- 500
chains <- 12

rerun_MCMC <- FALSE

Multivariate analysis

if (rerun_MCMC) {
  HS <- h2life %>% 
    filter(treat == "HS") %>% rename(Age_HS = NewAge) %>% 
    as.data.frame()
  DR <- h2life %>% 
    filter(treat == "DR") %>% rename(Age_DR = NewAge) %>% 
    as.data.frame()
  STD <- h2life %>%
    filter(treat == "STD") %>% rename(Age_STD = NewAge) %>% 
    as.data.frame()
  
  h2life_mrg <- full_join(full_join(HS, DR), STD)
  
  prior1 <- list(R = list(V = diag(3) * 1.002, nu = 1.002),
                 G = list(G1 = list(V = diag(3) * 1.002, nu = 0.002)))
  
  # prior2 <- list(R = list(V = diag(3) / 3, nu = 1.002),
  #                G = list(G1 = list(V = diag(3) / 3, nu = 0.002)))
  # 
  # prior3 <- list(R = list(V = diag(3) * 1.002, nu = 1),
  #                G = list(G1 = list(V = diag(3) * 0.5, nu = 0.002)))
  
  # priors <- list(prior1, prior2, prior3)
  # prior_names <- c("Tri: V = diag(3) * 1.002, nu = 0.02",
  #                  "Tri: V = diag(3) / 3, nu = 0.02",
  #                  "Tri: V = diag(3) * 0.5, nu = 0.02")
  priors <- list(prior1)
  prior_names <- c("Tri: V = diag(3) * 1.002, nu = 0.02")
  
  for (ii in 1:length(priors)) {
    prior <- priors[[ii]]
    fm <- mclapply(1:chains, function(i) {
      MCMCglmm(cbind(Age_STD, Age_HS, Age_DR) ~ trait - 1,
               random = ~ us(trait):animal,
               rcov = ~ idh(trait):units,
               data = h2life_mrg,
               prior = prior,
               pedigree = pedigree,
               family = rep("gaussian", 3),
               nitt = iter,
               burnin = burnin,
               thin = thinning)
    }, mc.cores = chains)
    outfile <- paste0("../../Data/Processed/surv_multivariate_model_prior", ii, ".Rda")
    save(fm, file = outfile)
    
    re <- lapply(fm, function(m) m$VCV)
    re <- do.call(mcmc.list, re)
    re <- as.mcmc(as.matrix(re))
    n_eff <- effectiveSize(re)
    
    # STD vs. HS
    HS_STD <- re[ , "traitAge_HS:traitAge_STD.animal"] /
      sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
             re[ , "traitAge_HS:traitAge_HS.animal"])
    
    # STD vs. DR
    DR_STD <- re[ , "traitAge_DR:traitAge_STD.animal"] /
      sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
             re[ , "traitAge_DR:traitAge_DR.animal"])
    
    # DR vs. HS
    HS_DR <- re[ , "traitAge_HS:traitAge_DR.animal"] /
      sqrt(re[ , "traitAge_DR:traitAge_DR.animal"] *
             re[ , "traitAge_HS:traitAge_HS.animal"])
    
    genet_corr <- bind_rows(genet_corr,
                            tibble(model = prior_names[[ii]],
                                   HS_STD = median(HS_STD),
                                   DR_STD = median(DR_STD),
                                   HS_DR = median(HS_DR),
                                   n_eff = mean(n_eff)))
  }
  
  write_csv(genet_corr, path = "../../Data/Processed/Lifespan_Genetic_Correlations.csv")
}

Analyze model

load("../../Data/Processed/surv_multivariate_model_prior1.Rda")

fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe[, 1, drop = FALSE], ask = FALSE)

plot(fe[, 2, drop = FALSE], ask = FALSE)

plot(fe[, 3, drop = FALSE], ask = FALSE)

# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)

plot(re[, 1, drop = FALSE], ask = FALSE)

plot(re[, 2, drop = FALSE], ask = FALSE)

plot(re[, 3, drop = FALSE], ask = FALSE)

autocorr.diag(re)
##           traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## Lag 0                          1.000000000                    1.0000000000
## Lag 500                        0.068179841                    0.0736391590
## Lag 2500                      -0.006882237                   -0.0038467280
## Lag 5000                       0.007044146                   -0.0061371156
## Lag 25000                      0.004727290                    0.0009293054
##           traitAge_DR:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## Lag 0                         1.000000000                    1.0000000000
## Lag 500                       0.063500507                    0.0736391590
## Lag 2500                     -0.004802190                   -0.0038467280
## Lag 5000                     -0.001234391                   -0.0061371156
## Lag 25000                    -0.001749820                    0.0009293054
##           traitAge_HS:traitAge_HS.animal traitAge_DR:traitAge_HS.animal
## Lag 0                        1.000000000                   1.0000000000
## Lag 500                      0.083782065                   0.0888369970
## Lag 2500                    -0.017718219                   0.0045442414
## Lag 5000                    -0.005514260                  -0.0107881960
## Lag 25000                   -0.001397025                  -0.0004649028
##           traitAge_STD:traitAge_DR.animal traitAge_HS:traitAge_DR.animal
## Lag 0                         1.000000000                   1.0000000000
## Lag 500                       0.063500507                   0.0888369970
## Lag 2500                     -0.004802190                   0.0045442414
## Lag 5000                     -0.001234391                  -0.0107881960
## Lag 25000                    -0.001749820                  -0.0004649028
##           traitAge_DR:traitAge_DR.animal traitAge_STD.units
## Lag 0                        1.000000000        1.000000000
## Lag 500                      0.085607455        0.059879186
## Lag 2500                     0.002057077       -0.007839362
## Lag 5000                     0.004328985       -0.001042706
## Lag 25000                    0.006550331        0.001383093
##           traitAge_HS.units traitAge_DR.units
## Lag 0          1.0000000000       1.000000000
## Lag 500        0.0665168077       0.062501759
## Lag 2500      -0.0165415030       0.008408074
## Lag 5000      -0.0014701485      -0.005139524
## Lag 25000     -0.0007576496       0.008627476
effectiveSize(re)
## traitAge_STD:traitAge_STD.animal  traitAge_HS:traitAge_STD.animal 
##                         21514.52                         20212.93 
##  traitAge_DR:traitAge_STD.animal  traitAge_STD:traitAge_HS.animal 
##                         22013.61                         20212.93 
##   traitAge_HS:traitAge_HS.animal   traitAge_DR:traitAge_HS.animal 
##                         19734.78                         19002.38 
##  traitAge_STD:traitAge_DR.animal   traitAge_HS:traitAge_DR.animal 
##                         22013.61                         19002.38 
##   traitAge_DR:traitAge_DR.animal               traitAge_STD.units 
##                         20088.25                         21512.98 
##                traitAge_HS.units                traitAge_DR.units 
##                         20320.53                         20642.51
# Concatenate samples
re <- as.mcmc(as.matrix(re))

head(re)
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1 
## End = 7 
## Thinning interval = 1 
##      traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## [1,]                        0.4632983                      0.16953957
## [2,]                        0.3950603                      0.13623892
## [3,]                        0.3263196                      0.01576263
## [4,]                        0.5508081                      0.12743270
## [5,]                        0.5017048                      0.12533453
## [6,]                        0.5403548                      0.25806069
## [7,]                        0.3794987                      0.10694896
##      traitAge_DR:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## [1,]                       0.1009902                      0.16953957
## [2,]                       0.1959820                      0.13623892
## [3,]                       0.1045252                      0.01576263
## [4,]                       0.1523802                      0.12743270
## [5,]                       0.1975044                      0.12533453
## [6,]                       0.1772361                      0.25806069
## [7,]                       0.1167529                      0.10694896
##      traitAge_HS:traitAge_HS.animal traitAge_DR:traitAge_HS.animal
## [1,]                      0.2790736                     0.11957617
## [2,]                      0.3379790                     0.10244792
## [3,]                      0.2690235                     0.07564472
## [4,]                      0.3205873                     0.06663020
## [5,]                      0.3030006                     0.07402833
## [6,]                      0.3937808                     0.13460634
## [7,]                      0.3468938                     0.12529748
##      traitAge_STD:traitAge_DR.animal traitAge_HS:traitAge_DR.animal
## [1,]                       0.1009902                     0.11957617
## [2,]                       0.1959820                     0.10244792
## [3,]                       0.1045252                     0.07564472
## [4,]                       0.1523802                     0.06663020
## [5,]                       0.1975044                     0.07402833
## [6,]                       0.1772361                     0.13460634
## [7,]                       0.1167529                     0.12529748
##      traitAge_DR:traitAge_DR.animal traitAge_STD.units traitAge_HS.units
## [1,]                      0.2258080          0.4949568         0.5178935
## [2,]                      0.2310788          0.5079044         0.4634037
## [3,]                      0.1762062          0.5622988         0.5508867
## [4,]                      0.1506864          0.4227779         0.4810127
## [5,]                      0.2671842          0.4244858         0.4844486
## [6,]                      0.1886153          0.4402332         0.4419436
## [7,]                      0.2341845          0.5260786         0.4923344
##      traitAge_DR.units
## [1,]         0.5405749
## [2,]         0.5475709
## [3,]         0.5945725
## [4,]         0.6132244
## [5,]         0.5262783
## [6,]         0.5440586
## [7,]         0.5453158
# STD vs. HS
plot(re[ , "traitAge_HS:traitAge_STD.animal"])

plot(re[ , "traitAge_STD:traitAge_STD.animal"])

HS_STD <- re[ , "traitAge_HS:traitAge_STD.animal"] /
  sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
         re[ , "traitAge_HS:traitAge_HS.animal"])
plot(HS_STD)

median(HS_STD)
## [1] 0.2999119
HPDinterval(HS_STD)
##           lower     upper
## var1 0.02156968 0.5526564
## attr(,"Probability")
## [1] 0.95
# STD vs. DR
DR_STD <- re[ , "traitAge_DR:traitAge_STD.animal"] /
  sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
         re[ , "traitAge_DR:traitAge_DR.animal"])
plot(DR_STD)

median(DR_STD)
## [1] 0.4957338
HPDinterval(DR_STD)
##          lower     upper
## var1 0.2584886 0.7007605
## attr(,"Probability")
## [1] 0.95
# DR vs. HS
HS_DR <- re[ , "traitAge_HS:traitAge_DR.animal"] /
  sqrt(re[ , "traitAge_DR:traitAge_DR.animal"] *
         re[ , "traitAge_HS:traitAge_HS.animal"])
plot(HS_DR)

median(HS_DR)
## [1] 0.3501346
HPDinterval(HS_DR)
##           lower     upper
## var1 0.07818288 0.6078474
## attr(,"Probability")
## [1] 0.95

Plot for paper

M <- data_frame(`HS vs. STD` = as.numeric(HS_STD),
                `DR vs. STD` = as.numeric(DR_STD),
                `HS vs. DR` = as.numeric(HS_DR))

M <- as_tibble( M ) %>%  
  select(`HS vs. STD`, `DR vs. STD`, `HS vs. DR`) %>% 
  rename(`HS vs. C` = `HS vs. STD`) %>% 
  rename(`DR vs. C` = `DR vs. STD`) %>% 
  rename(`HS vs. DR` = `HS vs. DR`)

colMeans(M)
##  HS vs. C  DR vs. C HS vs. DR 
## 0.2939628 0.4873072 0.3441595
M %>% gather(Comparison, value) %>%
  ggplot(aes(value, color = Comparison)) +
  geom_line(stat = "density", size = 1.0) +
  labs(x = "Genetic Correlation", y = "Density") +
  theme(legend.position = c(0.10, 0.85),
        text = element_text(size = 10),
        legend.text = element_text(size = 10)) +
  scale_x_continuous(limits = c(-1, 1)) +
  my_theme

ggsave(file = "../../Figures/Lifespan_Genetic_Correlation_Plot.pdf",
       width = 4, height = 4)
Lifespan_correlation <- last_plot()
save(Lifespan_correlation, file = "../../Figures/Lifespan_correlation.Rda")

Parameter expanded prior multivariate analysis

if (rerun_MCMC) {
  
  iter <- 1e5
  burnin <- 2e3
  thinning <- 500
  chains <- 12

  HS <- h2life %>% 
    filter(treat == "HS") %>% rename(Age_HS = NewAge) %>% 
    as.data.frame()
  DR <- h2life %>% 
    filter(treat == "DR") %>% rename(Age_DR = NewAge) %>% 
    as.data.frame()
  STD <- h2life %>%
    filter(treat == "STD") %>% rename(Age_STD = NewAge) %>% 
    as.data.frame()
  
  h2life_mrg <- full_join(full_join(HS, DR), STD)
  
  prior1 <- list(R = list(V = diag(3), nu = 1.002, fix = 3),
                 G = list(G1 = list(V = diag(3),
                                    nu = 1,
                                    alpha.mu = rep(0, 3),
                                    alpha.V = diag(3) * 25 ^ 2)))
  
  priors <- list(prior1)
  prior_names <- c("Tri: Parameter expanded")
  
  for (ii in 1:length(priors)) {
    prior <- priors[[ii]]
    fm <- mclapply(1:chains, function(i) {
      MCMCglmm(cbind(Age_STD, Age_HS, Age_DR) ~ trait - 1,
               random = ~ us(trait):animal,
               rcov = ~ idh(trait):units,
               data = h2life_mrg,
               prior = prior,
               pedigree = pedigree,
               family = rep("gaussian", 3),
               nitt = iter,
               burnin = burnin,
               thin = thinning)
    }, mc.cores = chains)
    outfile <- paste0("../../Data/Processed/surv_multivariate_model_prior_expanded", ii, ".Rda")
    save(fm, file = outfile)
    
    re <- lapply(fm, function(m) m$VCV)
    re <- do.call(mcmc.list, re)
    re <- as.mcmc(as.matrix(re))
    n_eff <- effectiveSize(re)
    
    # STD vs. HS
    HS_STD <- re[ , "traitAge_HS:traitAge_STD.animal"] /
      sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
             re[ , "traitAge_HS:traitAge_HS.animal"])
    
    # STD vs. DR
    DR_STD <- re[ , "traitAge_DR:traitAge_STD.animal"] /
      sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
             re[ , "traitAge_DR:traitAge_DR.animal"])
    
    # DR vs. HS
    HS_DR <- re[ , "traitAge_HS:traitAge_DR.animal"] /
      sqrt(re[ , "traitAge_DR:traitAge_DR.animal"] *
             re[ , "traitAge_HS:traitAge_HS.animal"])
    
    genet_corr <- bind_rows(genet_corr,
                            tibble(model = prior_names[[ii]],
                                   HS_STD = median(HS_STD),
                                   DR_STD = median(DR_STD),
                                   HS_DR = median(HS_DR),
                                   n_eff = mean(n_eff)))
  }
  
  write_csv(
    genet_corr,
    path = "../../Data/Processed/Lifespan_Genetic_Correlations_Parameter_Expanded.csv")
}

Analyze parameter expanded model

load("../../Data/Processed/surv_multivariate_model_prior_expanded1.Rda")

fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe[, 1, drop = FALSE], ask = FALSE)

plot(fe[, 2, drop = FALSE], ask = FALSE)

plot(fe[, 3, drop = FALSE], ask = FALSE)

# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)

plot(re[, 1, drop = FALSE], ask = FALSE)

plot(re[, 2, drop = FALSE], ask = FALSE)

plot(re[, 3, drop = FALSE], ask = FALSE)

autocorr.diag(re)
##           traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## Lag 0                          1.000000000                     1.000000000
## Lag 500                       -0.011880172                     0.084936923
## Lag 2500                      -0.012336107                    -0.044213068
## Lag 5000                      -0.008184613                     0.022182921
## Lag 25000                      0.009565219                     0.005296313
##           traitAge_DR:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## Lag 0                         1.000000000                     1.000000000
## Lag 500                       0.077587962                     0.084936923
## Lag 2500                     -0.010870186                    -0.044213068
## Lag 5000                      0.030051023                     0.022182921
## Lag 25000                    -0.007022395                     0.005296313
##           traitAge_HS:traitAge_HS.animal traitAge_DR:traitAge_HS.animal
## Lag 0                       1.0000000000                    1.000000000
## Lag 500                     0.0003537141                    0.014058983
## Lag 2500                    0.0136531239                   -0.023536549
## Lag 5000                   -0.0173873268                   -0.017871792
## Lag 25000                  -0.0131964670                    0.009574731
##           traitAge_STD:traitAge_DR.animal traitAge_HS:traitAge_DR.animal
## Lag 0                         1.000000000                    1.000000000
## Lag 500                       0.077587962                    0.014058983
## Lag 2500                     -0.010870186                   -0.023536549
## Lag 5000                      0.030051023                   -0.017871792
## Lag 25000                    -0.007022395                    0.009574731
##           traitAge_DR:traitAge_DR.animal traitAge_STD.units
## Lag 0                        1.000000000        1.000000000
## Lag 500                      0.003887223       -0.002186594
## Lag 2500                    -0.005213611       -0.024391116
## Lag 5000                    -0.033258672       -0.009846700
## Lag 25000                    0.001597150        0.006159378
##           traitAge_HS.units traitAge_DR.units
## Lag 0           1.000000000               NaN
## Lag 500         0.023515392               NaN
## Lag 2500       -0.001714462               NaN
## Lag 5000       -0.036410907               NaN
## Lag 25000      -0.008078434               NaN
effectiveSize(re)
## traitAge_STD:traitAge_STD.animal  traitAge_HS:traitAge_STD.animal 
##                         2339.831                         2209.375 
##  traitAge_DR:traitAge_STD.animal  traitAge_STD:traitAge_HS.animal 
##                         2331.718                         2209.375 
##   traitAge_HS:traitAge_HS.animal   traitAge_DR:traitAge_HS.animal 
##                         2737.938                         2323.729 
##  traitAge_STD:traitAge_DR.animal   traitAge_HS:traitAge_DR.animal 
##                         2331.718                         2323.729 
##   traitAge_DR:traitAge_DR.animal               traitAge_STD.units 
##                         2475.675                         2512.021 
##                traitAge_HS.units                traitAge_DR.units 
##                         2437.652                            0.000
# Concatenate samples
re <- as.mcmc(as.matrix(re))

head(re)
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1 
## End = 7 
## Thinning interval = 1 
##      traitAge_STD:traitAge_STD.animal traitAge_HS:traitAge_STD.animal
## [1,]                        0.3177375                    0.0001236148
## [2,]                        0.5171791                    0.0012672359
## [3,]                        0.3841529                    0.0197225670
## [4,]                        0.3423261                    0.0346445797
## [5,]                        0.3976341                    0.0676707225
## [6,]                        0.5314955                    0.1152669886
## [7,]                        0.3384343                    0.1062301494
##      traitAge_DR:traitAge_STD.animal traitAge_STD:traitAge_HS.animal
## [1,]                      0.08665067                    0.0001236148
## [2,]                      0.09078737                    0.0012672359
## [3,]                      0.07710740                    0.0197225670
## [4,]                      0.04821795                    0.0346445797
## [5,]                      0.06409478                    0.0676707225
## [6,]                      0.06619153                    0.1152669886
## [7,]                      0.05919508                    0.1062301494
##      traitAge_HS:traitAge_HS.animal traitAge_DR:traitAge_HS.animal
## [1,]                      0.2019285                     0.03867934
## [2,]                      0.3023297                     0.06524771
## [3,]                      0.2444371                     0.07475899
## [4,]                      0.2704081                     0.06776188
## [5,]                      0.2260605                     0.06126987
## [6,]                      0.2762569                     0.08347054
## [7,]                      0.2542423                     0.07237181
##      traitAge_STD:traitAge_DR.animal traitAge_HS:traitAge_DR.animal
## [1,]                      0.08665067                     0.03867934
## [2,]                      0.09078737                     0.06524771
## [3,]                      0.07710740                     0.07475899
## [4,]                      0.04821795                     0.06776188
## [5,]                      0.06409478                     0.06126987
## [6,]                      0.06619153                     0.08347054
## [7,]                      0.05919508                     0.07237181
##      traitAge_DR:traitAge_DR.animal traitAge_STD.units traitAge_HS.units
## [1,]                     0.07656928          0.5323402         0.5340417
## [2,]                     0.07191132          0.4571359         0.4873783
## [3,]                     0.08480175          0.4993428         0.5207448
## [4,]                     0.05734316          0.5560712         0.4897561
## [5,]                     0.07447087          0.5157063         0.5407287
## [6,]                     0.09450949          0.4034260         0.5127812
## [7,]                     0.10054775          0.5484401         0.5123482
##      traitAge_DR.units
## [1,]                 1
## [2,]                 1
## [3,]                 1
## [4,]                 1
## [5,]                 1
## [6,]                 1
## [7,]                 1
# STD vs. HS
plot(re[ , "traitAge_HS:traitAge_STD.animal"])

plot(re[ , "traitAge_STD:traitAge_STD.animal"])

HS_STD <- re[ , "traitAge_HS:traitAge_STD.animal"] /
  sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
         re[ , "traitAge_HS:traitAge_HS.animal"])
plot(HS_STD)

median(HS_STD)
## [1] 0.2311976
HPDinterval(HS_STD)
##            lower     upper
## var1 -0.03035282 0.5091342
## attr(,"Probability")
## [1] 0.9498299
# STD vs. DR
DR_STD <- re[ , "traitAge_DR:traitAge_STD.animal"] /
  sqrt(re[ , "traitAge_STD:traitAge_STD.animal"] *
         re[ , "traitAge_DR:traitAge_DR.animal"])
plot(DR_STD)

median(DR_STD)
## [1] 0.4468741
HPDinterval(DR_STD)
##          lower     upper
## var1 0.2074835 0.6598369
## attr(,"Probability")
## [1] 0.9498299
# DR vs. HS
HS_DR <- re[ , "traitAge_HS:traitAge_DR.animal"] /
  sqrt(re[ , "traitAge_DR:traitAge_DR.animal"] *
         re[ , "traitAge_HS:traitAge_HS.animal"])
plot(HS_DR)

median(HS_DR)
## [1] 0.3498687
HPDinterval(HS_DR)
##         lower     upper
## var1 0.112513 0.6110439
## attr(,"Probability")
## [1] 0.9498299

Plot parameter expanded model

M <- data_frame(`HS vs. STD` = as.numeric(HS_STD),
                `DR vs. STD` = as.numeric(DR_STD),
                `HS vs. DR` = as.numeric(HS_DR))

M <- as_tibble( M ) %>%  
  select(`HS vs. STD`, `DR vs. STD`, `HS vs. DR`) %>% 
  rename(`HS vs. C` = `HS vs. STD`) %>% 
  rename(`DR vs. C` = `DR vs. STD`) %>% 
  rename(`HS vs. DR` = `HS vs. DR`)

colMeans(M)
##  HS vs. C  DR vs. C HS vs. DR 
## 0.2304897 0.4403069 0.3474292
M %>% gather(Comparison, value) %>%
  ggplot(aes(value, color = Comparison)) +
  geom_line(stat = "density", size = 1.0) +
  labs(x = "Genetic Correlation", y = "Density") +
  theme(legend.position = c(0.10, 0.85),
        text = element_text(size = 10),
        legend.text = element_text(size = 10)) +
  scale_x_continuous(limits = c(-1, 1)) +
  my_theme

# ggsave(file = "../../Figures/Lifespan_Genetic_Correlation_Plot.pdf",
       # width = 4, height = 4)
# Lifespan_correlation <- last_plot()
# save(Lifespan_correlation, file = "../../Figures/Lifespan_correlation.Rda")